{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0a4a46ad",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using the Community license in this session. If you have a full Xpress license, first set the XPAUTH_PATH environment variable to the full path to your license file, xpauth.xpr, and then restart Python. If you want to use the FICO Community license and no longer want to see this message, set the XPAUTH_PATH environment variable to: /home/maged/anaconda3/envs/balance/lib/python3.8/site-packages/xpress/license/community-xpauth.xpr\n",
      "NB: setting XPAUTH_PATH will also affect any other Xpress products installed on your system.\n"
     ]
    }
   ],
   "source": [
    "# Importing modules\n",
    "try:\n",
    "    import hsbalance as hs;\n",
    "except ImportError:\n",
    "    !pip install hsbalance\n",
    "    import hsbalance as hs;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f622ade9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a602bc7e",
   "metadata": {},
   "source": [
    "# Huang test_case"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17f1f0c2",
   "metadata": {},
   "source": [
    "The configuration of the gas turbine we consider is shown in Fig. 1. This gas turbine has a mass of 727,000 Kg\n",
    "and is 61.65 m long.\n",
    "This gas turbine has three planes (BZ-A, BZ-C and BZ-E) available for balancing. However, in most cases, only\n",
    "two planes (BZ-A and BZ-E) are used because of the difficulty of placing weights in plane BZ-C. The objective\n",
    "of balancing here is to minimize the maximum vibration amplitude of the rotor at the measuring planes by using\n",
    "balance weights placed in planes BZ-A and BZ-E.\n",
    "For this case, the only available balance weight type for balance planes BZ-A and BZ-E is 142 grams. Each balance\n",
    "plane also has a limited number of weight holes, as shown in Fig. 1 (Plane BZ-A has balance holes placed every 7.5\n",
    "degrees, and plane BZ-E has balance holes placed every 5 degrees). It is possible to use all available balance holes.\n",
    "However, there is one extra constraint: at most one weight can be put in each balance hole.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f8088da8",
   "metadata": {},
   "outputs": [],
   "source": [
    "ALPHA_math=[\n",
    "                            ['0.085@27', '0.05@82'],\n",
    "                            ['0.053@57', '0.071@15']]\n",
    "\n",
    "A_math=[\n",
    "                            ['32@357'], \n",
    "                            ['105@346']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7144ede",
   "metadata": {},
   "source": [
    "Convert to complex numbers (cartesian) form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7933f461",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = hs.convert_matrix_to_cart(A_math)\n",
    "ALPHA = hs.convert_matrix_to_cart(ALPHA_math)\n",
    "# A, ALPHA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9524b2c9",
   "metadata": {},
   "source": [
    "Adding ALPHA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "eb13e858",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = hs.Alpha()\n",
    "alpha.add(direct_matrix=ALPHA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a41c59f0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/maged/Turbomachinery-Rotors-Balancing/src/hsbalance/IC_matrix.py:109: UserWarning: \n",
      "Warning: Influence Matrix is asymmetrical!\n",
      "  warnings.warn('\\nWarning: Influence Matrix is asymmetrical!')\n",
      "\n",
      "Influence Matrix is asymmetrical, check your data.\n",
      "\n",
      "No ill conditioned planes --> ok\n"
     ]
    }
   ],
   "source": [
    "alpha.check()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4fa0b512",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "Influence Coefficient Matrix\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "Coefficient Values\n",
      "==============================\n",
      "               Plane 1       Plane 2\n",
      "Sensor 1  0.085 @ 27.0   0.05 @ 82.0\n",
      "Sensor 2  0.053 @ 57.0  0.071 @ 15.0\n",
      "==============================\n",
      "End of Coefficient Values\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n"
     ]
    }
   ],
   "source": [
    "print(alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c81228e",
   "metadata": {},
   "source": [
    "## Solving with Least squares:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6b55e039",
   "metadata": {},
   "outputs": [],
   "source": [
    " # Instantiate least square model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2ec2321b",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_LeastSquares = hs.LeastSquares(A, alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "25d4e6ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "W_LeastSquares = model_LeastSquares.solve() #solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "918ed0cd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([['639.888 @ 73.8'],\n",
       "       ['1122.814 @ 165.2']], dtype='<U16')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hs.convert_matrix_to_math(W_LeastSquares)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "05b32e14",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_LeastSquares.rmse()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "663277a5",
   "metadata": {},
   "source": [
    "## Splitting plane BZ-A"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "668d510e",
   "metadata": {},
   "source": [
    "* Holes available: we select angles from 30-150 to reduce calculation time, holes spaced by 7.5 degrees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ed534718",
   "metadata": {},
   "outputs": [],
   "source": [
    "split_BZA = model_LeastSquares.create_split()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "3fe906a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "angles = np.arange(30,150,7.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6edfacb8",
   "metadata": {},
   "source": [
    "* Weights available is 142 gram (5 ounces) and only one weight per hole is allowed.\n",
    "* We can also constrain the problem to 2000 grams per plane to avoid solution infiltration. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "596fbc86",
   "metadata": {},
   "outputs": [],
   "source": [
    "split_BZA.split_setup(0, max_number_weights_per_hole=1, holes_available=[angles]\n",
    "                         ,weights_available=[142], max_weight_per_plane=2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bb1d377",
   "metadata": {},
   "source": [
    "Benchmarking the split_solve method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "2dc3155c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "50.2 ms ± 3.22 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "split_BZA.split_solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d1956af",
   "metadata": {},
   "source": [
    "Solving the split problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "ea1b1cab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 1., 0., 1., 0., 0., 1., 0., 0., 1., 1., 0., 0., 0., 0., 0.]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZA.split_solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e48f5f4",
   "metadata": {},
   "source": [
    "Getting the results in pandas dataframe (including none zero values only)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "9559b8d7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>37.5</th>\n",
       "      <th>52.5</th>\n",
       "      <th>75.0</th>\n",
       "      <th>97.5</th>\n",
       "      <th>105.0</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>weights_available</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>142</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  37.5  52.5  75.0  97.5  105.0\n",
       "weights_available                              \n",
       "142                 1.0   1.0   1.0   1.0   1.0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZA.results()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "589ea6b1",
   "metadata": {},
   "source": [
    "df returns dataframe that describes the number of weights(index) for every hole (columns).  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9df4d213",
   "metadata": {},
   "source": [
    "Calculating the error weight (unbalance weight value) caused by splitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "bff59a8b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.10034846])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZA.error()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b83a1363",
   "metadata": {},
   "source": [
    "Calculating the equivalent weight caused by splitting "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "32a45124",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array('640.2 @ 73.6', dtype='<U12')"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZA.error(options='equ')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a143cfd6",
   "metadata": {},
   "source": [
    "which is close to the weight generated by the model '639.888@73.8'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "49f03026",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZA.error(options='problem_error') # getting the problem accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4473353b",
   "metadata": {},
   "source": [
    "  - prob_error: As I tried to solve the problem using only free solvers, I used [`mip_cvxpy`][mip_cvxpy] package that uses `CBC` solver to solver mip problems. To get the most accurate results we should have turned the problem into mixed integer second order cone which is not supported by CBC, so instead I reformulate the objective function to be linear optimizing the max or real and imaginary part of the result weight error. The splitting error is returned then so the user can estimate.\n",
    "  - I finally settled on using solver = `XPRESS` which is free with limited matrix value of 5000. I believe in practice this is satisfying for balancing problem.\n",
    "      - [`FICO® Xpress`][xpress] solver is [fast and accurate](./test/benchmark_splitting_solver.md) as it allows quadratic solving for the MIP problem.\n",
    "[xpress]: https://www.fico.com/en/products/fico-xpress-solver\n",
    "[mip_cvxpy]: https://pypi.org/project/mip-cvxpy/"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24470248",
   "metadata": {},
   "source": [
    "We can update the original W_leastSquares weights with the new split weight we have just created"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "0994402c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/maged/Turbomachinery-Rotors-Balancing/src/hsbalance/model.py:256: UserWarning: This will change your model optimum solution.Choose confirm=True\n",
      "  warnings.warn('This will change your model optimum solution.'\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[  178.48435221+614.4911797j ],\n",
       "       [-1085.53108741+286.93648042j]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZA.update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "48da341b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([['639.888 @ 73.8'],\n",
       "       ['1122.814 @ 165.2']], dtype='<U16')"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hs.convert_matrix_to_math(W_LeastSquares)  # The new W_LeastSquares"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f074d84",
   "metadata": {},
   "source": [
    "By default update will not do anything but warning you, as this will change the optimum solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "1966f443",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  180.56557795+614.20840222j],\n",
       "       [-1085.53108741+286.93648042j]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZA.update(confirm=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "51399f94",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([['640.2 @ 73.6'],\n",
       "       ['1122.814 @ 165.2']], dtype='<U16')"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hs.convert_matrix_to_math(W_LeastSquares)  # The new W_LeastSquares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "01cfa4bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1449"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_LeastSquares.rmse()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae39818a",
   "metadata": {},
   "source": [
    "We can see the RMS error has increased from 0 to 0.145 in this problem due to splitting. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a458177e",
   "metadata": {},
   "source": [
    "Print a report of the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "47129eef",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "MODEL\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "MODEL TYPE\n",
      "==================================================\n",
      "LeastSquares\n",
      "==================================================\n",
      "End of MODEL TYPE\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "INFLUENCE COEFFICIENT MATRIX\n",
      "==================================================\n",
      "\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "Influence Coefficient Matrix\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "Coefficient Values\n",
      "==============================\n",
      "               Plane 1       Plane 2\n",
      "Sensor 1  0.085 @ 27.0   0.05 @ 82.0\n",
      "Sensor 2  0.053 @ 57.0  0.071 @ 15.0\n",
      "==============================\n",
      "End of Coefficient Values\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "==================================================\n",
      "End of INFLUENCE COEFFICIENT MATRIX\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "INITIAL VIBRATION\n",
      "==================================================\n",
      "              Vibration\n",
      "Sensor 1   32.0 @ 357.0\n",
      "Sensor 2  105.0 @ 346.0\n",
      "==================================================\n",
      "End of INITIAL VIBRATION\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "SOLUTION\n",
      "==================================================\n",
      "        Correction Masses\n",
      "Plane 1      640.2 @ 73.6\n",
      "Plane 2  1122.814 @ 165.2\n",
      "==================================================\n",
      "End of SOLUTION\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "RMSE\n",
      "==================================================\n",
      "0.1449\n",
      "==================================================\n",
      "End of RMSE\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "Resiudal Vibration Expected\n",
      "==================================================\n",
      "         Expected Vibration\n",
      "Sensor 1       0.179 @ 19.3\n",
      "Sensor 2       0.111 @ 49.3\n",
      "==================================================\n",
      "End of Resiudal Vibration Expected\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "SPLITS\n",
      "==================================================\n",
      "                  37.5  52.5  75.0  97.5  105.0\n",
      "weights_available                              \n",
      "142                 1.0   1.0   1.0   1.0   1.0\n",
      "==================================================\n",
      "End of SPLITS\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n"
     ]
    }
   ],
   "source": [
    "print(model_LeastSquares.info())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ef457ce",
   "metadata": {},
   "source": [
    "## Splitting plane BZ-C"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71e39a7d",
   "metadata": {},
   "source": [
    "* Holes available: we select angles from 100-200 to reduce calculation time, holes spaced by 5 degrees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "7dc70508",
   "metadata": {},
   "outputs": [],
   "source": [
    "split_BZC = model_LeastSquares.create_split()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "3da7b0dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "angles = np.arange(100, 200, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ae13706",
   "metadata": {},
   "source": [
    "* Weights available is 142 gram (5 ounces) and only one weight per hole is allowed.\n",
    "* We can also constrain the problem to 2000 grams per plane to avoid solution infiltration. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "656c2e54",
   "metadata": {},
   "outputs": [],
   "source": [
    "split_BZC.split_setup(1, max_number_weights_per_hole=1, holes_available=[angles]\n",
    "                         ,weights_available=[142], max_weight_per_plane=2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a310499a",
   "metadata": {},
   "source": [
    "Benchmarking the split_solve method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "e147a0cf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "93.7 ms ± 9.19 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "split_BZC.split_solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43df6795",
   "metadata": {},
   "source": [
    "Solving the split problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "4cf6b9f2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0., 0., 0., 0., 0., 1., 0., 0., 1., 0., 0., 1., 0., 1.,\n",
       "        1., 1., 1., 1.]])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZC.split_solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "323370e6",
   "metadata": {},
   "source": [
    "Getting the results in pandas dataframe (including none zero values only)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "f3a6ee65",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr:last-of-type th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>100</th>\n",
       "      <th>135</th>\n",
       "      <th>150</th>\n",
       "      <th>165</th>\n",
       "      <th>175</th>\n",
       "      <th>180</th>\n",
       "      <th>185</th>\n",
       "      <th>190</th>\n",
       "      <th>195</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>weights_available</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>142</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   100  135  150  165  175  180  185  190  195\n",
       "weights_available                                             \n",
       "142                1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZC.results()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64beba01",
   "metadata": {},
   "source": [
    "df returns dataframe that describes the number of weights(index) for every hole (columns).  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1d29c31",
   "metadata": {},
   "source": [
    "Calculating the error weight (unbalance weight value) caused by splitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "20a4b752",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.63300893])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZC.error()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df4d0ec4",
   "metadata": {},
   "source": [
    "Calculating the equivalent weight caused by splitting "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "e4afe820",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array('1124.27 @ 165.2', dtype='<U15')"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZC.error(options='equ')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cd803d6",
   "metadata": {},
   "source": [
    "which is close to the weight generated by the model '1122.814@165.2'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "f2e72abd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.87276186e-16])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZC.error(options='problem_error') # getting the problem accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3941c957",
   "metadata": {},
   "source": [
    "  - prob_error: As I tried to solve the problem using only free solvers, I used [`mip_cvxpy`][mip_cvxpy] package that uses `CBC` solver to solver mip problems. To get the most accurate results we should have turned the problem into mixed integer second order cone which is not supported by CBC, so instead I reformulate the objective function to be linear optimizing the max or real and imaginary part of the result weight error. The splitting error is returned then so the user can estimate.\n",
    "  - I finally settled on using solver = `XPRESS` which is free with limited matrix value of 5000. I believe in practice this is satisfying for balancing problem.\n",
    "      - [`FICO® Xpress`][xpress] solver is [fast and accurate](./test/benchmark_splitting_solver.md) as it allows quadratic solving for the MIP problem.\n",
    "[xpress]: https://www.fico.com/en/products/fico-xpress-solver\n",
    "[mip_cvxpy]: https://pypi.org/project/mip-cvxpy/"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c2bd06b",
   "metadata": {},
   "source": [
    "We can update the original W_leastSquares weights with the new split weight we have just created"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "64bf90e5",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/maged/Turbomachinery-Rotors-Balancing/src/hsbalance/model.py:256: UserWarning: This will change your model optimum solution.Choose confirm=True\n",
      "  warnings.warn('This will change your model optimum solution.'\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[  180.56557795+614.20840222j],\n",
       "       [-1085.53108741+286.93648042j]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZC.update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "baf89158",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([['640.2 @ 73.6'],\n",
       "       ['1122.814 @ 165.2']], dtype='<U16')"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hs.convert_matrix_to_math(W_LeastSquares)  # The new W_LeastSquares"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2334f63b",
   "metadata": {},
   "source": [
    "By default update will not do anything but warn you, as this will change the optimum solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "779ee7b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  180.56557795+614.20840222j],\n",
       "       [-1087.12774135+286.59382263j]])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "split_BZC.update(confirm=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "11ae469a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([['640.2 @ 73.6'],\n",
       "       ['1124.27 @ 165.2']], dtype='<U15')"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hs.convert_matrix_to_math(W_LeastSquares)  # The new W_LeastSquares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "0e6d9308",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1099"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_LeastSquares.rmse()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad6eff97",
   "metadata": {},
   "source": [
    "Print a report of the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "5374e943",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "MODEL\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "MODEL TYPE\n",
      "==================================================\n",
      "LeastSquares\n",
      "==================================================\n",
      "End of MODEL TYPE\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "INFLUENCE COEFFICIENT MATRIX\n",
      "==================================================\n",
      "\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "Influence Coefficient Matrix\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "Coefficient Values\n",
      "==============================\n",
      "               Plane 1       Plane 2\n",
      "Sensor 1  0.085 @ 27.0   0.05 @ 82.0\n",
      "Sensor 2  0.053 @ 57.0  0.071 @ 15.0\n",
      "==============================\n",
      "End of Coefficient Values\n",
      "++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "==================================================\n",
      "End of INFLUENCE COEFFICIENT MATRIX\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "INITIAL VIBRATION\n",
      "==================================================\n",
      "              Vibration\n",
      "Sensor 1   32.0 @ 357.0\n",
      "Sensor 2  105.0 @ 346.0\n",
      "==================================================\n",
      "End of INITIAL VIBRATION\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "SOLUTION\n",
      "==================================================\n",
      "        Correction Masses\n",
      "Plane 1      640.2 @ 73.6\n",
      "Plane 2   1124.27 @ 165.2\n",
      "==================================================\n",
      "End of SOLUTION\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "RMSE\n",
      "==================================================\n",
      "0.1099\n",
      "==================================================\n",
      "End of RMSE\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "Resiudal Vibration Expected\n",
      "==================================================\n",
      "         Expected Vibration\n",
      "Sensor 1      0.176 @ 352.6\n",
      "Sensor 2      0.044 @ 134.1\n",
      "==================================================\n",
      "End of Resiudal Vibration Expected\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "SPLITS\n",
      "==================================================\n",
      "                  37.5  52.5  75.0  97.5  105.0\n",
      "weights_available                              \n",
      "142                 1.0   1.0   1.0   1.0   1.0\n",
      "\n",
      "                   100  135  150  165  175  180  185  190  195\n",
      "weights_available                                             \n",
      "142                1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0\n",
      "==================================================\n",
      "End of SPLITS\n",
      "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n",
      "\n",
      "                   \n"
     ]
    }
   ],
   "source": [
    "print(model_LeastSquares.info())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84511dd7",
   "metadata": {},
   "source": [
    "We can see the rms error has increased from 0 to 0.1099 in this problem due to splitting. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "2773f6ad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<hsbalance.model._Model.Split at 0x7fece465d970>,\n",
       " <hsbalance.model._Model.Split at 0x7fed78097e50>]"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_LeastSquares.split_instance   # We can see the splitting instances that override the optimum solution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03102a20",
   "metadata": {},
   "source": [
    "### Discussion\n",
    "- The results in the Huang paper are different from mine. But I would like to further prove that the above calculations give even much lesser error.  \n",
    "- The time is an issue in splitting, you should solve the problem first and try to select suitable angles around the original phase angle."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ad4b4e0",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baabbb2e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7410adcc",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:balance] *",
   "language": "python",
   "name": "conda-env-balance-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
